In [88]:
from IPython.display import HTML

HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Mostra codice/ Nascondi codice."></form>''')
Out[88]:
In [89]:
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
from sklearn import metrics
import pydotplus
import matplotlib.colors as mcolors
import statsmodels.formula.api as smf 
from sklearn.metrics import mean_squared_error
from IPython.core.display import HTML
from sklearn.model_selection import train_test_split, cross_val_score, cross_val_predict, StratifiedKFold
import plotly.express as px
from sklearn.model_selection import train_test_split, GridSearchCV, cross_val_score, RepeatedStratifiedKFold, StratifiedKFold
from sklearn.metrics import accuracy_score, confusion_matrix,roc_curve, roc_auc_score, precision_score, recall_score, precision_recall_curve
from sklearn.metrics import f1_score

import warnings
warnings.filterwarnings('ignore')
from sklearn.metrics import fbeta_score
from sklearn.metrics import f1_score
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score
from sklearn.svm import SVC
from sklearn.preprocessing import OneHotEncoder
from sklearn.model_selection import StratifiedShuffleSplit
from sklearn import preprocessing  # to normalise existing X
from sklearn.cluster import KMeans
from sklearn.preprocessing import LabelEncoder
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn import preprocessing as prp
import plotly.express as px
import ipywidgets as widgets
import plotly.graph_objs as go
from sklearn.linear_model import LinearRegression
import statsmodels.formula.api as stat
from sklearn.model_selection import train_test_split
lm =LinearRegression()
from plotly.offline import plot, iplot
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import classification_report
from sklearn.metrics import confusion_matrix,accuracy_score

Homework II- Business Analytics- Simona Maria Borrello - S277789

In [3]:
data = pd.read_csv('C:/Users/Simona/OneDrive/Desktop/MAGISTRALE/Secondo semestre/Business Analytc/casi/pilgA/PilgrimABC.csv', sep=',')

Controlliamo le dimensioni del dataset e verifichiamo la presenza di valori mancanti nei dati.

In [4]:
print('osservazioni, variabili =', data.shape[0], 'x', data.shape[1])
print('Dati mancanti per variabile')
print(data.isnull().sum())
osservazioni, variabili = 31634 x 11
Dati mancanti per variabile
ID              0
9Profit         0
9Online         0
9Age         8289
9Inc         8261
9Tenure         0
9District       0
0Profit      5238
0Online      5219
9Billpay        0
0Billpay     5219
dtype: int64

Decidiamo di eliminare le colonne '0Billpay', '9Billpay' per la nostra analisi e introduciamo alcune nuove variabili al fine di gestire i valori mancanti per gli attributi '9Age', '9Inc', '0Profit' e '0Online'.

In particolare introduciamo le seguenti nuove variabili:

  • 'AgeGiven', che vale 1 se l'età del cliente è nota;
  • 'Retain', che vale 1 se il cliente ha abbandonato la banca;
  • ' AgeZero', che assume lo stesso valore di '9Age' per i dati presenti e 0 per quelli mancanti;
  • 'IncomeGiven', che vale 1 se è noto il reddito del cliente:
  • 'Profit0', che assegna 0 ai dati mancanti di 'Profit0';
  • 'IncomeZero', che assegna 0 ai dati mancanti in '9Inc';
  • 'AgeAverage' , che assume lo stesso valore di '9Age' per i dati presenti e la media degli altri valori per quelli mancanti.
In [5]:
data = data.drop(columns=['0Billpay', '9Billpay'])

data['AgeGiven']=  np.where(data['9Age'].isnull(), 0, 1) #1 ho l'informazione della età
data['Retain']=  np.where(data['0Profit'].isnull(), 1, 0) #0 ho conservato il cliente
data['AgeZero'] = data['9Age'].fillna(0)
data['IncomeGiven']=  np.where(data['9Inc'].isnull(), 0, 1) #1 ho l'informazione della data
data['Profit0'] = data['0Profit'].fillna(0)
data['IncomeZero'] = data['9Inc'].fillna(0)
data['AgeAverage'] = data['9Age'].fillna(data['9Age'].mean()) 
data['Online0'] = data['0Online'].fillna(0)

data.rename(columns={"9Profit": "Profit9", "9Online": "Online9", "9Tenure": "Tenure9", "9District":"District9"}, inplace=True)

#data.describe()
data2= data.copy()
data2= data2.drop(['ID','9Age','9Inc','0Online', '0Profit'], axis=1)
data=data2.copy()

Ecco il dataset, dopo le modifiche:

In [6]:
data
Out[6]:
Profit9 Online9 Tenure9 District9 AgeGiven Retain AgeZero IncomeGiven Profit0 IncomeZero AgeAverage Online0
0 21 0 6.33 1200 0 1 0.0 0 0.0 0.0 4.046048 0.0
1 -6 0 29.50 1200 1 0 6.0 1 -32.0 3.0 6.000000 0.0
2 -49 1 26.41 1100 1 0 5.0 1 -22.0 5.0 5.000000 1.0
3 -4 0 2.25 1200 0 1 0.0 0 0.0 0.0 4.046048 0.0
4 -61 0 9.91 1200 1 0 2.0 1 -4.0 9.0 2.000000 0.0
... ... ... ... ... ... ... ... ... ... ... ... ...
31629 -50 0 3.75 1200 1 0 5.0 1 1.0 5.0 5.000000 0.0
31630 458 0 12.08 1300 1 0 3.0 1 423.0 8.0 3.000000 1.0
31631 -83 0 15.83 1200 1 0 6.0 1 -60.0 4.0 6.000000 0.0
31632 92 1 5.41 1200 1 0 1.0 1 170.0 6.0 1.000000 1.0
31633 124 0 17.50 1300 1 0 3.0 1 150.0 6.0 3.000000 0.0

31634 rows × 12 columns

Codifichiamo le variabili categoriche e verifichiamo che il dataset sia bilanciato, rispetto alla variabile 'Retain'.

In [8]:
columns_to_encode = ['District9','AgeZero','IncomeZero']
encoder=OneHotEncoder(sparse=False)
data_pre= data

X_encoded = pd.DataFrame (encoder.fit_transform(data_pre[['District9']]))
X_encoded.columns = encoder.get_feature_names(['District9'])
data_pre.drop(['District9'] ,axis=1, inplace=True)
p1= pd.concat([data_pre, X_encoded ], axis=1)

X_encoded = pd.DataFrame (encoder.fit_transform(p1[['AgeZero']]))
X_encoded.columns = encoder.get_feature_names(['AgeZero'])
p1.drop(['AgeZero'] ,axis=1, inplace=True)
p2= pd.concat([p1, X_encoded ], axis=1)


X_encoded = pd.DataFrame (encoder.fit_transform(p2[['IncomeZero']]))
X_encoded.columns = encoder.get_feature_names(['IncomeZero'])
p2.drop(['IncomeZero'] ,axis=1, inplace=True)
data_enc= pd.concat([p2, X_encoded ], axis=1)

data_enc.rename(columns={'AgeZero_1.0':'AgeZero_1fascia',
                          'AgeZero_2.0':'AgeZero_2fascia',
                          'AgeZero_3.0':'AgeZero_3fascia',
                           'AgeZero_0.0':'AgeZero_0fascia',
                            'AgeZero_4.0':'AgeZero_4fascia',
                             'AgeZero_5.0':'AgeZero_5fascia',
                              'AgeZero_6.0':'AgeZero_6fascia',
                              'AgeZero_7.0':'AgeZero_7fascia'}, 
                 inplace=True)

data_enc.rename(columns={'IncomeZero_1.0':'IncomeZero_1fascia',
                          'IncomeZero_2.0':'IncomeZero_2fascia',
                          'IncomeZero_3.0':'IncomeZero_3fascia',
                           'IncomeZero_0.0':'IncomeZero_0fascia',
                            'IncomeZero_4.0':'IncomeZero_4fascia',
                             'IncomeZero_5.0':'IncomeZero_5fascia',
                              'IncomeZero_6.0':'IncomeZero_6fascia',
                              'IncomeZero_7.0':'IncomeZero_7fascia', 
                              'IncomeZero_8.0':'IncomeZero_8fascia',
                              'IncomeZero_9.0':'IncomeZero_9fascia'}, 
                 inplace=True)

data_enc.rename(columns={'District9_1300.0':'District_thir',
                          'District9_1200.0':'District_sec',
                          
                              'District9_1100.0':'District_first'}, 
                 inplace=True)
In [9]:
data_clas= data_enc.drop(columns=['Online0'])
In [10]:
group_size=data['Retain'].value_counts().sort_index()
group_names=['Clienti Fedeli (' + '{:1.1f}'.format(group_size[0]/(group_size[0]+group_size[1])*100) + '%)',
             'Clienti Persi (' + '{:1.1f}'.format(group_size[1]/(group_size[0]+group_size[1])*100) + '%)']

a, b=[plt.cm.Greens, plt.cm.Reds]

fig, ax = plt.subplots()
ax.axis('equal')

mypie, _ = ax.pie(group_size,  labels=group_names, colors=[a(0.6), b(0.6)], explode=(0.1,0),
                  textprops={'size': 18, 'color':'black'})
plt.figsize=(5,5)

plt.setp( mypie, width=1, edgecolor='white')
plt.margins(0,0)
fig = plt.gcf()
fig.set_size_inches(4,4) 

plt.show()

Come si nota il dataset è sbilanciato, questo implica la necessità di utilizzare, in fase di classificazione, determinate tecniche che tengano conto della distribuzione del dataset.

Visualizziamo il dataset mediante alcuni grafici, ad esempio la distribuzione delle diverse fasce d'età nel presenti nel dataset.

In [11]:
a=31634
#vediamo se è bilanciato

group_size=data['AgeAverage'].value_counts().sort_index()
group_names=['Fascia-età1 (' + '{:1.1f}'.format(group_size[1]/(a)*100) + '%)',
             'Fascia-età2 (' + '{:1.1f}'.format(group_size[2]/(a)*100) + '%)',
             'Fascia-età3 (' + '{:1.1f}'.format(group_size[3]/(a)*100) + '%)', 
             'Fascia-età4 (' + '{:1.1f}'.format(group_size[4]/(a)*100) + '%)', 
             'Età assente (' + '{:1.1f}'.format(8289/(a)*100) + '%)',
             'Fascia-età5 (' + '{:1.1f}'.format(group_size[5]/(a)*100) + '%)',
             'Fascia-età6 (' + '{:1.1f}'.format(group_size[6]/(a)*100) + '%)',
             'Fascia-età7 (' + '{:1.1f}'.format(group_size[7]/(a)*100) + '%)']
#colorlist = ['crimson','k','r','gold','g','b','c','d','a','l']
colors = mcolors.TABLEAU_COLORS.values()

fig, ax = plt.subplots()

ax.axis('equal')

mypie, _ = ax.pie( group_size, radius=2.5, labels=group_names, colors=colors,
                  textprops={'size': 18, 'color':'black'})
plt.setp( mypie, width=0.7, edgecolor='white')
plt.margins(0,0)
fig.set_size_inches(3,3) 

plt.show()

#grafic
In [12]:
#grafico2
x1 = data_enc.loc[data_enc.AgeGiven==0, 'Profit9']
x2 = data_enc.loc[data_enc.AgeGiven==1, 'Profit9']

fig, axes = plt.subplots(1, 2, figsize=(10, 3), sharey=True, dpi=100)
axes[0].title.set_text('Dato su età mancante')
axes[1].title.set_text('Dato su età presente')

sns.distplot(x1 , color="dodgerblue", ax=axes[0], axlabel='Profitto9')
sns.distplot(x2 , color="deeppink", ax=axes[1], axlabel='Profitto9')
#plt.hist(x1, alpha=0.5, bins=100, density=True, stacked=True, label='age', color='dodgerblue')
#plt.hist(x2, alpha=0.5, bins=100, density=True,stacked=True, label='noage', color='deeppink')


plt.tight_layout();
In [13]:
barWidth = 0.25
plt.figure(figsize=(10,5)) 


# set width of bar
barWidth = 0.4
a= data_enc.loc[data_enc.IncomeGiven==0, 'AgeAverage']
b= data_enc.loc[data_enc.IncomeGiven==1, 'AgeAverage']
a= group_size=a.value_counts().sort_index()
b= group_size=b.value_counts().sort_index()

# set height of bar
bars1 = a
bars2 = b
 
# Set position of bar on X axis
r1 = np.arange(len(bars1))
r2 = [x + barWidth for x in r1]
 
# Make the plot
plt.bar(r1, bars1, color='g', width=barWidth, edgecolor='white', label=' INCOME mancante')
plt.bar(r2, bars2, color='orange', width=barWidth, edgecolor='white', label=' INCOME presente')

# Add xticks on the middle of the group bars
plt.xticks([r + barWidth for r in range(len(bars1))],  ['Age1', 'Age2', 'Age3', 'Age4', 'Età mancante','Age5','Age6','Age7'], fontweight='bold', fontsize=20)
plt.xticks(rotation=90)
plt.yticks(fontsize=20)
# Create legend & Show graphic


plt.legend(fontsize=15)
plt.show()

Il grafico precedente mostra che nella maggior parte delle osservazioni in cui non è presente il dato su Income, non è presente neppure il dato su Età.

In [14]:
sns.set(rc={'figure.figsize':(15,10)})
sns.set(font_scale = 2)

bplot=sns.boxplot(y='Profit9', x='Retain', 
                 data=data_enc, 
                 width=0.5,
                 palette="nipy_spectral")

Il boxplot mostra, in entrambe le classi ( Retain=1 o Retain=0) , la presenza di numerosi (eventuali) outliers .

In [15]:
data= data2.copy()

Regressione

Cerchiamo di capire, grazie ad alcune semplici regressioni lineari, in che modo le variabili 'Age' , 'Income' e 'Online'influenzino il profitto .

Ruolo della variabile 'Age'

Decidiamo come trattare la variabile Age. Abbiamo diverse alternative:

  • come variabile categorica, a cui associamo la fascia 0 ai dati in cui l'età è mancante
  • come variabile numerica:
    • in cui la media è sostituita ai dati mancanti ;
    • in cui 0 è sostituito ai dati mancanti .
In [16]:
X= data.drop(['Retain', 'Profit9','Profit0'], axis=1)
Y= data['Profit9']
data_reg= X.join(Y)

Proviamo a vedere la relazione tra 'Profit9' e 'AgeGiven'.

In [22]:
modello1 = smf.ols('Profit9 ~ AgeGiven', data=data_enc).fit()
In [23]:
def short_summary(est):
    return HTML(est.summary().tables[1].as_html())

print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello1.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello1.rsquared_adj))
short_summary(modello1)
Il valore di MSE del modello vale 73916.26, mentre l'R^2adj vale 0.0071
Out[23]:
coef std err t P>|t| [0.025 0.975]
Intercept 72.9625 2.986 24.433 0.000 67.109 78.816
AgeGiven 52.2245 3.476 15.024 0.000 45.411 59.038

Il modello precedente mostra che la variabile 'AgeGiven è significativa (p-value < 0.05) . In particolare il coefficiente è positivo per cui ci aspettiamo che il profitto per i dati di cui è nota l'età sia superiore, in media , rispetto al profitto delle osservazioni di cui non è nota l'età.

Infatti se vogliamo entrare più in dettaglio, il valore di intercetta non è altro che la media della variabile 'Profit9' per i dati in cui 'AgeGiven' è uguale a 0.

In [24]:
a= data_enc[data_enc['AgeGiven']==0]
a["Profit9"].mean()
Out[24]:
72.96248039570516

Invece, il coefficiente $\beta$ di AgeGiven non è altro che la differenza tra l'intercetta e la media della variabile 'Profit9' per i dati in cui 'AgeGiven' è uguale a 1.

In [25]:
b= data_enc[data_enc['AgeGiven']==1]
b["Profit9"].mean()- a["Profit9"].mean()
Out[25]:
52.224497543896476

Adesso, proviamo a vedere la relazione tra 'Profit9' e la variabile 'Age' come categorica.

In [26]:
modello1 = smf.ols('Profit9 ~  AgeZero_1fascia+ AgeZero_2fascia+ AgeZero_3fascia+ AgeZero_4fascia+ AgeZero_5fascia+ AgeZero_6fascia + AgeZero_7fascia', data=data_enc).fit()
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello1.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello1.rsquared_adj))
short_summary(modello1)
#rint(results.summary())
#Y_pred= modello1.predict(X)
#mse = mean_squared_error(Y, Y_pred)


#_pred= results.predict(X_test)
#mse = mean_squared_error(Y_test, Y_pred)

#results = smf.ols('Profit9 ~  Online9+  AgeGiven', data=data_reg).fit()
#print(results.summary())
Il valore di MSE del modello vale 72532.01, mentre l'R^2adj vale 0.0256
Out[26]:
coef std err t P>|t| [0.025 0.975]
Intercept 72.9625 2.958 24.665 0.000 67.164 78.760
AgeZero_1fascia -69.5076 10.531 -6.600 0.000 -90.149 -48.866
AgeZero_2fascia -14.4729 5.350 -2.705 0.007 -24.959 -3.987
AgeZero_3fascia 42.1498 4.712 8.944 0.000 32.913 51.386
AgeZero_4fascia 62.6993 4.716 13.295 0.000 53.455 71.943
AgeZero_5fascia 72.7971 5.583 13.040 0.000 61.855 83.739
AgeZero_6fascia 87.4476 6.358 13.754 0.000 74.986 99.909
AgeZero_7fascia 119.2989 5.974 19.971 0.000 107.590 131.007

Come si nota il valore dell'intercetta di questo modello è uguale al valore dell'intercetta del modello precedente. Questo non ci stupisce, perchè in questo caso la categoria di riferimento della regressione è proprio la classe di clienti di cui non è nota l'età. Dai coefficienti della regressione abbiamo una informazione in più sulle medie della variabile 'Profit9'. Infatti i coefficienti sono negativi per le categorie 1 e 2, il chè implica che la media delle variabile 'Profit9' è minore rispetto alla categoria di riferimento.

Mentre i clienti delle restanti categorie sono, in media, più profittevoli rispetto a coloro di cui non è nota l'età. A titolo esemplificativo, vediamo il calcolo del coefficiente relativo a AgeZero_1fascia:

In [27]:
b= data_enc[data_enc['AgeZero_1fascia']==1]
b["Profit9"].mean()- a["Profit9"].mean()
Out[27]:
-69.50755081824036

Adesso proviamo a considerare la variabile 'AgeAverage' ( numerica) .

In [28]:
modello2 = smf.ols('Profit9 ~  AgeAverage', data=data_reg).fit()
In [29]:
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello2.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello2.rsquared_adj))
short_summary(modello2)
Il valore di MSE del modello vale 73205.65, mentre l'R^2adj vale 0.0166
Out[29]:
coef std err t P>|t| [0.025 0.975]
Intercept 10.5967 4.620 2.293 0.022 1.541 19.653
AgeAverage 24.9394 1.078 23.129 0.000 22.826 27.053

Infine consideriamo la variabile 'AgeZero' (numerica) .

In [30]:
modello3 = smf.ols('Profit9 ~  AgeZero', data=data_reg).fit()
In [31]:
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello3.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello3.rsquared_adj))
short_summary(modello3)
Il valore di MSE del modello vale 72855.14, mentre l'R^2adj vale 0.0213
Out[31]:
coef std err t P>|t| [0.025 0.975]
Intercept 59.0939 2.507 23.571 0.000 54.180 64.008
AgeZero 17.5523 0.668 26.262 0.000 16.242 18.862

Scegliamo di trattare la variabile sull'età come categorica sia perchè sembra essere la scelta più ragionevole, confermata anche dal valore di MSE e R^2adj dei modelli precedenti.

Ruolo della variabile 'Income'

Vediamo la relazione tra 'Profit9' e 'IncomeGiven'.

In [32]:
modello1 = smf.ols('Profit9 ~ IncomeGiven', data=data_enc).fit()
In [33]:
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello1.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello1.rsquared_adj))
short_summary(modello1)
Il valore di MSE del modello vale 73874.24, mentre l'R^2adj vale 0.0076
Out[33]:
coef std err t P>|t| [0.025 0.975]
Intercept 71.3650 2.990 23.865 0.000 65.504 77.226
IncomeGiven 54.3241 3.479 15.615 0.000 47.505 61.143

La variabile 'IncomeGiven' sembra essere significativa ed, in particolare, il profitto è, in media , maggiore per le osservazioni in cui è presente il dato su 'Income'.

Adesso, proviamo a vedere la relazione tra 'Profit9' e la variabile 'Income' come categorica.

In [34]:
modello1 = smf.ols('Profit9 ~  IncomeZero_1fascia+ IncomeZero_2fascia+ IncomeZero_3fascia+ IncomeZero_4fascia + IncomeZero_5fascia+ IncomeZero_6fascia+ IncomeZero_7fascia+ IncomeZero_8fascia+ IncomeZero_9fascia', data=data_enc).fit()
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello1.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello1.rsquared_adj))
short_summary(modello1)
Il valore di MSE del modello vale 72243.82, mentre l'R^2adj vale 0.0295
Out[34]:
coef std err t P>|t| [0.025 0.975]
Intercept 71.3650 2.957 24.132 0.000 65.569 77.161
IncomeZero_1fascia 0.9129 6.640 0.137 0.891 -12.102 13.928
IncomeZero_2fascia 18.0881 9.896 1.828 0.068 -1.309 37.485
IncomeZero_3fascia 17.5942 6.070 2.899 0.004 5.697 29.492
IncomeZero_4fascia 22.0148 6.324 3.481 0.000 9.620 34.410
IncomeZero_5fascia 23.6540 6.264 3.776 0.000 11.376 35.932
IncomeZero_6fascia 47.9035 4.700 10.192 0.000 38.691 57.116
IncomeZero_7fascia 66.6487 5.627 11.844 0.000 55.619 77.678
IncomeZero_8fascia 86.4628 7.086 12.201 0.000 72.573 100.352
IncomeZero_9fascia 162.5114 5.758 28.225 0.000 151.226 173.797

Guardando i coefficienti notiamo che essi sono crescenti, fatta eccezione per la terza fascia. Si potrebbe pensare che il valore di Profitto è, in qualche modo, proporzionale al valore di Income.

Ruolo della variabile 'Online'

Cerchiamo di capire, in che modo, il fatto che il cliente utilizzi o meno il servizio Online, incida sul profitto.

In [35]:
#Grafico
on= data[data['Online9']==1]
off= data[data['Online9']==0]
plt.figure(figsize=(15,5))
ax = sns.regplot(x="Profit9", y="Profit0", data=on)
plt.title('Utenti che usano i servizi Online', fontsize=25)
plt.xlabel('Profit_1999', fontsize=20)
plt.ylabel('Profit_2000', fontsize=2)
plt.tick_params( labelsize=15)
In [36]:
plt.figure(figsize=(15,5))
#plt.plot(off['Profit9'], off['Profit0'], marker="o", color='green', linestyle="")


ax = sns.regplot(x="Profit9", y="Profit0", data=off, color='green')
#Grafico
plt.title('Utenti che non usano i servizi Online', fontsize=25)
plt.xlabel('Profit_1999', fontsize=20)
plt.ylabel('Profit_2000', fontsize=20)
plt.tick_params( labelsize=15)

Vediamo la relazione tra 'Profit9' e 'Online9'.

In [38]:
modello2 = smf.ols('Profit9 ~  Online9 ', data=data_reg).fit()
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello2.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello2.rsquared_adj))
short_summary(modello2)
Il valore di MSE del modello vale 74439.99, mentre l'R^2adj vale 0.0000
Out[38]:
coef std err t P>|t| [0.025 0.975]
Intercept 110.7862 1.637 67.678 0.000 107.578 113.995
Online9 5.8806 4.690 1.254 0.210 -3.312 15.073

Vediamo la relazione tra 'Profit0' e 'Online0' ( per i clienti che non hanno abbandonato la banca).

In [39]:
modello2 = smf.ols('Profit0 ~  Online0', data=data).fit()
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello2.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello2.rsquared_adj))
short_summary(modello2)
Il valore di MSE del modello vale 129492.41, mentre l'R^2adj vale 0.0024
Out[39]:
coef std err t P>|t| [0.025 0.975]
Intercept 112.8539 2.216 50.935 0.000 108.511 117.197
Online0 48.1039 5.436 8.850 0.000 37.450 58.758

Per entrambi gli anni, il profitto è, in media, maggiore per i clienti che usano il servizio online.

Inoltre cerchiamo di capire come cambia, in media, il profitto, se un cliente passa da Online9=0 a Online0=1, cioè per i clienti che iniziano a usufrire del servizio online. Per tale analisi, considereremo, per ovvie ragione, solo i clienti che non hanno abbandonato la banca.

L'obiettivo è capire se la variazione di profitto di un cliente sia attribuibile, in parte, al fatto che il cliente abbia iniziato ad usare il servizio Online. Inoltre introduciamo una nuova colonna, che non è altro che la differenza tra i valori di Profit0 e Profit9.

Definiamo anche una nuova variabile ('start') come $ Online0 \cdot (1- Online9) $. Il valore di questa variabile assume 1 solo per i clienti che iniziano ad usare il servizio Online al secondo anno. Infatti :

  • Online9=0 , Online0=0 --> start= 0
  • Online9=0 , Online0=1 --> start =1
  • Online9=1 , Online0=0--> start =0
  • Online9=1 , Online0=1 --> start =0

Riassumendo i clienti con 'start'=1 sono i clienti che nel 1999 non usano il servizio Online e che nel 2000 lo usano.

In [40]:
data_online= data[data['Retain']==0]
data_online['diff'] = data_online['Profit0'] - data_online['Profit9']
data_online['start'] = data_online['Online0'] *(1- data_online['Online9'])
data_online
Out[40]:
Profit9 Online9 Tenure9 District9 AgeGiven Retain AgeZero IncomeGiven Profit0 IncomeZero AgeAverage Online0 diff start
1 -6 0 29.50 1200 1 0 6.0 1 -32.0 3.0 6.000000 0.0 -26.0 0.0
2 -49 1 26.41 1100 1 0 5.0 1 -22.0 5.0 5.000000 1.0 27.0 0.0
4 -61 0 9.91 1200 1 0 2.0 1 -4.0 9.0 2.000000 0.0 57.0 0.0
5 -38 0 2.33 1300 0 0 0.0 1 14.0 3.0 4.046048 0.0 52.0 0.0
6 -19 0 8.41 1300 1 0 3.0 1 0.0 1.0 3.000000 0.0 19.0 0.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
31629 -50 0 3.75 1200 1 0 5.0 1 1.0 5.0 5.000000 0.0 51.0 0.0
31630 458 0 12.08 1300 1 0 3.0 1 423.0 8.0 3.000000 1.0 -35.0 1.0
31631 -83 0 15.83 1200 1 0 6.0 1 -60.0 4.0 6.000000 0.0 23.0 0.0
31632 92 1 5.41 1200 1 0 1.0 1 170.0 6.0 1.000000 1.0 78.0 0.0
31633 124 0 17.50 1300 1 0 3.0 1 150.0 6.0 3.000000 0.0 26.0 0.0

26396 rows × 14 columns

Per i dati che soddisfano ' start= 1' la media della differenza tra i profitti vale:

In [41]:
a= data_online[data_online['start']==1]
a['diff'].mean()
Out[41]:
44.380554333479985
In [42]:
plt.figure(figsize=(13,7))
plt.plot(a['Profit9'], a['Profit0'], marker="o", color='brown', linestyle="")

plt.title('Utenti che hanno iniziato ad  usare i servizi Online', fontsize=25)
plt.xlabel('Profit_1999', fontsize=20)
plt.ylabel('Profit_2000', fontsize=20)
plt.tick_params( labelsize=15)

Per i dati che soddisfano ' start= 0' la media della differenza tra i profitti vale:

In [43]:
a= data_online[data_online['start']==0]
a['diff'].mean()
Out[43]:
23.180533101189734
In [44]:
plt.figure(figsize=(13,7))
plt.plot(a['Profit9'], a['Profit0'], marker="o", color='green', linestyle="")

plt.title('Utenti con start=0', fontsize=25)
plt.xlabel('Profit_1999', fontsize=20)
plt.ylabel('Profit_2000', fontsize=20)
plt.tick_params( labelsize=15)

Vogliamo stabilire se la differenza delle due medie è staticamente significativa. Per farlo, utlizziamo il seguente modello: $$ diff= Profit_0- Profit_9 = \beta_0 + \beta \cdot start $$

In [46]:
modello2 = smf.ols('diff ~ start', data=data_online).fit()
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello2.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello2.rsquared_adj))
short_summary(modello2)
Il valore di MSE del modello vale 99473.96, mentre l'R^2adj vale 0.0003
Out[46]:
coef std err t P>|t| [0.025 0.975]
Intercept 23.1805 2.031 11.415 0.000 19.200 27.161
start 21.2000 6.920 3.064 0.002 7.636 34.764

Vediamo che l'intercetta è la media della differenza tra i profitti dei due anni per i clienti con 'start=0' ( che abbiamo trovato precedentemente). Poichè il p-value del coefficiente 'start' vale 0.002, possiamo affermare, che la differenza tra le due medie è significativa e cioè rifiutiamo l'ipotesi che il coefficiente di 'start' sia zero.

Inoltre, dato che il coefficiente di 'start' è positivo , si potrebbe pensare che risulta conveniente per la banca che il cliente inizi ad usare il servizio online.

Verifichiamo cosa succede per i clienti che hanno smesso di usare il servizio online ( e che non hanno abbandonato la banca). Introduciamo la variabile 'stop'= $ Online9 \cdot (1- Online0) $. Il valore di questa variabile assume 1 solo per i clienti che hanno smesso di usare il servizio Online nel 2000. Infatti :

  • Online9=0 , Online0=0 --> stop= 0
  • Online9=0 , Online0=1 --> stop = 0
  • Online9=1 , Online0=0--> stop = 1
  • Online9=1 , Online0=1 --> stop =0

Quello che mi aspetterei è vedere il coefficiente della regressione negativo.

In [47]:
data_online['stop'] = data_online['Online9'] *(1- data_online['Online0'])

Calcoliamo la media della differenza tra i profitti per i clienti che hanno smesso di usare il servizio Online.

In [48]:
a= data_online[data_online['stop']==1]
a['diff'].mean()
Out[48]:
47.00308641975309

Calcoliamo la media della differenza tra i profitti per i clienti con 'Stop=0'.

In [49]:
a= data_online[data_online['stop']==0]
a['diff'].mean()
Out[49]:
24.73274010432648

Contro le mie aspettative, la media della differenza tra i profitti è maggiore per coloro che hanno abbandonato il servizio Online. Ma vediamo se la differenza tra le medie è statisticamente significativa.

In [51]:
modello2 = smf.ols('diff ~ stop', data=data_online).fit()
print('Il valore di MSE del modello vale' + ' {:1.2f}'.format(modello2.mse_resid)+ ', mentre l\'R^2adj vale' + ' {:1.4f}'.format(modello2.rsquared_adj))
short_summary(modello2)
Il valore di MSE del modello vale 99503.32, mentre l'R^2adj vale 0.0000
Out[51]:
coef std err t P>|t| [0.025 0.975]
Intercept 24.7327 1.954 12.660 0.000 20.904 28.562
stop 22.2703 17.633 1.263 0.207 -12.291 56.832

Quello che vediamo è che il coefficiente della regressione è positivo, questo indica che la media per i clienti che hanno smesso di usare il servizio online è più alta ( come abbiamo verificato sopra).

Tuttavia, in questo caso, il p-value del coefficiente 'stop' vale 0.207 > 0.05 e non possiamo rifiutare l'ipotesi che il coefficiente sia zero..

In [52]:
b= data_online[data_online['stop']==1]

plt.figure(figsize=(13,7))
plt.plot(b['Profit9'], b['Profit0'], marker="o", color='brown', linestyle="")

plt.title('Utenti che hanno smesso di usare i servizi Online', fontsize=25)
plt.xlabel('Profit_1999', fontsize=20)
plt.ylabel('Profit_2000', fontsize=20)
plt.tick_params( labelsize=15)

Regressione di 'Profit9'

Eliminazione Outliers

Avendo a disposizione un grande numero di dati, per trovare gli eventuali outliers, calcoliamo gli Z-score delle variabili numeriche e rimuoviamo le osservazioni che hanno un valore di Z-score, in valore assoluto, maggiore di 4.

In particolare, dato che per la regressione di 'Profit9' non utilizziamo come regressore la variabile 'Profit0' ( lo scopo è predire il valore di Profit9 quindi assumiamo di non conoscere il valore di profitto di un anno successivo), le uniche variabili numeriche di cui calcoliamo lo Z-score sono 'Profit9' e 'Tenure9'. Escludiamo dai regressori per la previsione di Profit9 anche la variabile 'Online0'.

Vediamo, mediante i seguenti istrogrammi, quante sono le variabili con uno Z-score maggiore, in valore assoluto, di 4.

In [53]:
import plotly.graph_objects as go

from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
#ciao=data_enc[data_enc['Retain']==0]
ciao= data_enc.copy()
l= ciao.shape
salva= ciao.copy()
scaled = scaler.fit_transform(ciao[['Profit9','Tenure9', 'Profit0']])
scaled = pd.DataFrame(scaled, columns=['Profit9','Tenure9', 'Profit0'])


plt.title('Profit9', fontsize=25)

plt.hist(scaled['Profit9'], bins = 50)
plt.xlabel('Zscore', fontweight='bold', fontsize=20)
plt.xticks(np.arange(-2, max(scaled['Profit9'])+1, 1.0), fontsize=20)
plt.yticks(fontsize=20)

plt.figure(figsize=(20,10))


plt.title('Tenure9', fontsize=25)

plt.hist(scaled['Tenure9'], bins = 50)
plt.xlabel('Zscore', fontweight='bold', fontsize=20)
plt.xticks(np.arange(-2, max(scaled['Tenure9'])+1, 1.0), fontsize=20 )
plt.yticks(fontsize=20)
plt.figure(figsize=(20,10))
Out[53]:
<Figure size 1440x720 with 0 Axes>
<Figure size 1440x720 with 0 Axes>
In [55]:
tol=4
ciao= data_enc.copy()
d = np.logical_or(abs(scaled['Tenure9'])>tol,abs(scaled['Profit9'])>tol)
ciao=pd.concat((ciao,d.rename('col')), axis=1, join='inner')
ciao=ciao[ciao['col']==False]
data_enc=ciao.drop(columns=['col'])
c=data_enc.shape
print("Il numero di osservazioni eliminate è",  (l[0]- c[0]))
Il numero di osservazioni eliminate è 384

Dividiamo il dataset in dataset train e test (70%- 30%).

Forward selection

Proviamo a predire la variabile 'Profit9'.

In [66]:
data_enc
X= data_enc.drop(['Profit9','Retain', 'Profit0','AgeAverage','Online0'], axis=1)
Y= data_enc['Profit9']
#Y= np.power(Y,2)
X, X_test, Y, Y_test = train_test_split(X, Y, test_size=0.3, random_state=42)

data_reg= X.join(Y)

ciao=[]
n = X.columns.shape[0] 

X_train, X_val, Y_train, Y_val = train_test_split(X, Y, test_size=0.5, random_state=42)

selected_features = [] 
left_feats = X.columns.tolist() 

for k in range(1,n+1):
    max_r_sq = 0.0
    chosen_feat = ""

    for feat in left_feats: 
        model = LinearRegression()
        model.fit(X_train[selected_features+[feat]],Y_train)
        r_sq = model.score(X_train[selected_features+[feat]],Y_train)
        
        if r_sq > max_r_sq: #if a better model is found (according to R2)
            chosen_feat = feat
            max_r_sq = r_sq
    
    selected_features.append(chosen_feat) 
    left_feats.remove(chosen_feat)       
print('Le variabili sono selezionate ( in base al valore di R^2) nel seguente ordine:\n')
j=1
for feat in selected_features:
    print( '-' + feat)
    j +=1
print('\n\n')

best_mse = -50000000000000
best_k = 0
for k in range(1,n+1):
    mse = cross_val_score(LinearRegression(), X_val[selected_features[:k]],
                          Y_val, cv=4, scoring='neg_mean_squared_error')
    mse = sum(mse)/5.0 #finding average model mse

    if mse>best_mse:
        best_mse = mse
        best_k = k

print('\n\nIl miglior modello, scelto sulla base del valore di MSE della cross validation (K=5) , è quello con ' + str(best_k) + ' variabili e i coefficienti stimati sono :\n')
model = LinearRegression().fit(X[selected_features[:best_k]],Y)
coefs = model.coef_
Y_pred = model.predict(X[selected_features[:best_k]])

res = Y-Y_pred
q= Y_pred
Y_preds = model.predict(X_test[selected_features[:best_k]])
mse = mean_squared_error(Y_test, Y_preds)
resA= res
yA= Y

for j in range(best_k):
    print(selected_features[j].ljust(20,' ') + '{:10.3f}'.format(coefs[j]))
print('\nIl TEST MSE del modello vale ' + '{:1.3f}'.format(mse) + '.')
Le variabili sono selezionate ( in base al valore di R^2) nel seguente ordine:

-Tenure9
-IncomeZero_9fascia
-AgeZero_7fascia
-IncomeZero_8fascia
-IncomeZero_7fascia
-IncomeZero_6fascia
-AgeZero_1fascia
-AgeZero_2fascia
-District_sec
-IncomeZero_1fascia
-AgeZero_6fascia
-Online9
-AgeZero_4fascia
-AgeZero_3fascia
-IncomeGiven
-AgeGiven
-IncomeZero_2fascia
-IncomeZero_5fascia
-District_thir
-IncomeZero_3fascia
-District_first
-AgeZero_0fascia
-IncomeZero_0fascia
-AgeZero_5fascia
-IncomeZero_4fascia





Il miglior modello, scelto sulla base del valore di MSE della cross validation (K=5) , è quello con 18 variabili e i coefficienti stimati sono :

Tenure9                  3.901
IncomeZero_9fascia     105.487
AgeZero_7fascia         54.863
IncomeZero_8fascia      68.808
IncomeZero_7fascia      47.984
IncomeZero_6fascia      27.794
AgeZero_1fascia        -56.850
AgeZero_2fascia        -25.096
District_sec            16.221
IncomeZero_1fascia     -19.387
AgeZero_6fascia         12.037
Online9                 15.536
AgeZero_4fascia         -0.790
AgeZero_3fascia          3.096
IncomeGiven             -2.081
AgeGiven                 1.807
IncomeZero_2fascia     -12.369
IncomeZero_5fascia       5.761

Il TEST MSE del modello vale 48769.177.

Regressione di Profit0

Eliminazione Outliers

Come dataset per le seguenti regressioni, consideriamo solo le osservazioni in cui la variabile Retain vale 0, cioè se il cliente non ha abbandonato la banca. Avendo già elimininato gli outliers rispetto alla variabile 'Profit9' e 'Tenure9', non ci resta che eliminarle rispetto a 'Profit0'.

In [58]:
data_enc= data_enc[data_enc['Retain']==0]
data_reg= X.join(Y)
data_enc
l=data_enc.shape

scaler = StandardScaler()
ciao= data_enc.copy()
salva= ciao.copy()
scaled = scaler.fit_transform(ciao[['Profit0']])
scaled = pd.DataFrame(scaled, columns=['Profit0'])
colors = px.colors.qualitative.Plotly





plt.title('Profit0', fontsize=25)

plt.hist(scaled['Profit0'], bins =50)
plt.xlabel('Zscore', fontweight='bold', fontsize=20)
plt.xticks(np.arange(-4, max(scaled['Profit0'])+1, 10) )
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.figure(figsize=(10,7))
Out[58]:
<Figure size 720x504 with 0 Axes>
<Figure size 720x504 with 0 Axes>
In [59]:
tol=4
ciao= data_enc.copy()

a= abs(scaled['Profit0'])>tol
ciao=pd.concat((ciao,a.rename('col')), axis=1, join='inner')
ciao=ciao[ciao['col']==False]
data_enc=ciao.drop(columns=['col'])
c=data_enc.shape
#print("Il numero di osservazioni eliminate è",  (l[0]- c[0]))

#c[0]

Dopo aver eliminato gli outliers, nuovamente, dividiamo il dataset in train e test.

In [60]:
X= data_enc.drop(['Profit9','Retain', 'Profit0','AgeAverage'], axis=1)
Y= data_enc['Profit0']

X, X_test, Y, Y_test = train_test_split(X, Y, test_size=0.3, random_state=42)

data_reg= X.join(Y)

Forward selection

Esclusione di 'Profit9' dai regressori

Proviamo a prevedere Profit0 senza usare Profit9.

In [61]:
ciao=[]
n = X.columns.shape[0] #max number of variables in the model


#cv
#sampling training set and validation set
X_train, X_val, Y_train, Y_val = train_test_split(X, Y, test_size=0.5, random_state=42)

selected_features = [] #empty list where to put selected features
left_feats = X.columns.tolist() #still to be chosen features

for k in range(1,n+1):
    max_r_sq = 0.0
    chosen_feat = ""

    for feat in left_feats: #looking for the next feature to select
        model = LinearRegression()
        model.fit(X_train[selected_features+[feat]],Y_train)
        r_sq = model.score(X_train[selected_features+[feat]],Y_train)
        
        if r_sq > max_r_sq: #if a better model is found (according to R2)
            chosen_feat = feat
            max_r_sq = r_sq
    
    selected_features.append(chosen_feat) #adding best feature to selected features list
    left_feats.remove(chosen_feat)        #and removing it from those still to be chosen
print('Le variabili sono selezionate ( in base al valore di R^2) nel seguente ordine:\n')
j=1
for feat in selected_features:
    print( '-' + feat)
    j +=1
print('\n\n')

#cross-validating and finding the best model
best_mse = -50000000000000
best_k = 0
for k in range(1,n+1):
    mse = cross_val_score(LinearRegression(), X_val[selected_features[:k]],
                          Y_val, cv=4, scoring='neg_mean_squared_error')
    mse = sum(mse)/5.0 #finding average model mse

    if mse>best_mse:
        best_mse = mse
        best_k = k

print('\n\nIl miglior modello, scelto sulla base del valore di MSE della cross validation ( K=5) , è quello con  ' + str(best_k) + ' variabili e i coefficienti stimati sono:\n')
model = LinearRegression().fit(X[selected_features[:best_k]],Y)
coefs = model.coef_
Y_pred = model.predict(X[selected_features[:best_k]])
#mse = mean_squared_error(Y, Y_pred)

res = Y-Y_pred
Y_preds = model.predict(X_test[selected_features[:best_k]])
mse = mean_squared_error(Y_test, Y_preds)


for j in range(best_k):
    print(selected_features[j].ljust(20,' ') + '{:10.3f}'.format(coefs[j]))
print('\nIl  TEST MSE del modello vale ' + '{:1.3f}'.format(mse) + '.')
Le variabili sono selezionate ( in base al valore di R^2) nel seguente ordine:

-Tenure9
-IncomeZero_9fascia
-IncomeZero_8fascia
-AgeZero_7fascia
-IncomeZero_1fascia
-Online0
-AgeZero_3fascia
-AgeGiven
-IncomeZero_7fascia
-IncomeZero_6fascia
-AgeZero_1fascia
-AgeZero_2fascia
-IncomeZero_3fascia
-AgeZero_6fascia
-District_first
-IncomeZero_2fascia
-AgeZero_4fascia
-Online9
-District_sec
-IncomeGiven
-IncomeZero_4fascia
-District_thir
-AgeZero_0fascia
-AgeZero_5fascia
-IncomeZero_0fascia
-IncomeZero_5fascia





Il miglior modello, scelto sulla base del valore di MSE della cross validation ( K=5) , è quello con  18 variabili e i coefficienti stimati sono:

Tenure9                  4.360
IncomeZero_9fascia     128.500
IncomeZero_8fascia     100.310
AgeZero_7fascia         82.283
IncomeZero_1fascia     -23.460
Online0                 24.623
AgeZero_3fascia         29.294
AgeGiven               -34.801
IncomeZero_7fascia      63.609
IncomeZero_6fascia      28.979
AgeZero_1fascia        -38.948
AgeZero_2fascia         -8.536
IncomeZero_3fascia       3.501
AgeZero_6fascia         29.164
District_first         -21.176
IncomeZero_2fascia     -31.903
AgeZero_4fascia         16.179
Online9                  8.781

Il  TEST MSE del modello vale 82047.749.

Inclusione di 'Profit9' nei regressori

Proviamo a predire profit0 usando profit9

In [62]:
X= data_enc.drop(['Retain', 'Profit0', 'Online0' ,'AgeAverage'], axis=1)
Y= data_enc['Profit0']

X, X_test, Y, Y_test = train_test_split(X, Y, test_size=0.3, random_state=42)

data_reg= X.join(Y)

ciao=[]
n = X.columns.shape[0] #max number of variables in the model


#cv
#sampling training set and validation set
X_train, X_val, Y_train, Y_val = train_test_split(X, Y, test_size=0.5, random_state=42)

selected_features = [] #empty list where to put selected features
left_feats = X.columns.tolist() #still to be chosen features

for k in range(1,n+1):
    max_r_sq = 0.0
    chosen_feat = ""

    for feat in left_feats: #looking for the next feature to select
        model = LinearRegression()
        model.fit(X_train[selected_features+[feat]],Y_train)
        r_sq = model.score(X_train[selected_features+[feat]],Y_train)
        
        if r_sq > max_r_sq: #if a better model is found (according to R2)
            chosen_feat = feat
            max_r_sq = r_sq
    
    selected_features.append(chosen_feat) #adding best feature to selected features list
    left_feats.remove(chosen_feat)        #and removing it from those still to be chosen
print('Le variabili sono selezionate ( in base al valore di R^2) nel seguente ordine:\n')
j=1
for feat in selected_features:
    print( '-' + feat)
    j +=1
print('\n\n')

#cross-validating and finding the best model
best_mse = -50000000000000
best_k = 0
for k in range(1,n+1):
    mse = cross_val_score(LinearRegression(), X_val[selected_features[:k]],
                          Y_val, cv=4, scoring='neg_mean_squared_error')
    mse = sum(mse)/5.0 #finding average model mse

    if mse>best_mse:
        best_mse = mse
        best_k = k
print('\n\nIl miglior modello, scelto sulla base del valore di MSE della cross validation ( K=5) , è quello con  ' + str(best_k) + ' variabili e i coefficienti stimati sono:\n')

model = LinearRegression().fit(X[selected_features[:best_k]],Y)
coefs = model.coef_
Y_pred = model.predict(X[selected_features[:best_k]])
msetr = mean_squared_error(Y, Y_pred)

res = Y-Y_pred

Y_preds = model.predict(X_test[selected_features[:best_k]])
mse = mean_squared_error(Y_test, Y_preds)


for j in range(best_k):
    print(selected_features[j].ljust(20,' ') + '{:10.3f}'.format(coefs[j]))
print('\nIl TEST MSE del modello  è ' + '{:1.3f}'.format(mse) + '.')
Le variabili sono selezionate ( in base al valore di R^2) nel seguente ordine:

-Profit9
-IncomeZero_8fascia
-IncomeZero_9fascia
-AgeZero_5fascia
-AgeGiven
-AgeZero_3fascia
-AgeZero_7fascia
-Online9
-IncomeZero_6fascia
-District_thir
-IncomeZero_7fascia
-IncomeZero_3fascia
-IncomeZero_2fascia
-AgeZero_4fascia
-Tenure9
-IncomeGiven
-AgeZero_2fascia
-IncomeZero_4fascia
-AgeZero_1fascia
-District_first
-IncomeZero_1fascia
-District_sec
-AgeZero_0fascia
-AgeZero_6fascia
-IncomeZero_0fascia
-IncomeZero_5fascia





Il miglior modello, scelto sulla base del valore di MSE della cross validation ( K=5) , è quello con  4 variabili e i coefficienti stimati sono:

Profit9                  0.866
IncomeZero_8fascia      36.159
IncomeZero_9fascia      33.838
AgeZero_5fascia        -20.681

Il TEST MSE del modello  è 49227.630.

Osservazioni

  • La presenza di Profit9 come regressore per la previsione di 'Profit0' fa quasi dimezzare il TEST MSE , determinando anche una notevole diminuizione del numero di regressori ( da 11 a 4). Questo fatto non stupisce perchè è chiaro che conoscore il valore del Profitto di un anno è determinante nella previsione del profitto dell'anno successivo. Notiamo, come ci si può aspettare, che la prima variabile selezionata dalla feature selection sia proprio Profit9.

  • La maggior parte delle variabili selezionate per la previsione 'Profit9' viene usata anche per la previsione 'Profit0' ( primo caso) . Anche questo fatto non stupisce in quanto ci si può aspettare che i fattori determinanti nella previsione del profitto non cambino al trascorrere degli anni.

Analisi dei residui

Riportiamo l'analisi dei residui per il primo dei modelli sopra analizzati.

In [63]:
plt.figure(figsize=(10,6))

ax= sns.scatterplot(x=yA, y=resA) 
ax.set_xlabel('Profit9 ')
ax.set_ylabel('Residui')

plt.show()
plt.figure(figsize=(10,8))

p= sns.scatterplot(x=q, y=resA) 
p.set_xlabel('Valori predetti ( fitted values )')
p.set_ylabel('Residui')

plt.show()

I grafici mostrano che le ipotesi alla base della regressione lineare non sono rispettate. Infatti è evidente un pattern nel grafico dei residui. Ho provato a effettuare alcune trasformazioni delle variabili. In particolare, ho provato ad elevare a diverse potenze le variabili numeriche del dataset, mentre ho escluso la trasformazione logaritmica perchè le variabili assumono anche valore negativo. Queste trasformazioni non hanno tuttavia determinato un miglioramento dei residui.

Di seguito, vediamo i grafici dei residui per la previsione del quadrato di Profit9.

In [65]:
plt.figure(figsize=(10,6))

ax= sns.scatterplot(x=yA, y=resA) 
ax.set_xlabel('(Profit9)^2')
ax.set_ylabel('Residui')

plt.show()
plt.figure(figsize=(10,8))

p= sns.scatterplot(x=q, y=resA) 
p.set_xlabel('Valori predetti ( fitted values )')
p.set_ylabel('Residui')

plt.show()

Classificazione

L'obiettivo è quello di assegnare ai dati una specifica classe ( positiva o negativa) , in base al valore di Retain. Ricordiamo che Retain vale:

  • 1 ,se il cliente ha abbandonato la banca -> etichetta positivo ( evento anomalo);
  • 0 , se il cliente resta in banca -> etichetta negativo.

In particolare, per la nostra analisi, assumiamo che i due possibili errori di classificazione abbiano un peso diverso. Consideriamo, infatti, che l'errore 'falso negativo' (il cliente abbandona la banca , ma viene classificato come un cliente che resta in banca) abbia una maggiore importanza rispetto a 'falso positivo'. Per sintetizzare, l'obiettivo dei modelli che seguono sarà quello di minimizzare i falsi negativi.

Inoltre, per prima cosa dividiamo il dataset in 'train data' e 'test data'. Quest'ultima porzione di dati verrà utilizzata esclusivamente alla fine della nostra analisi come strumento di confronto tra i diversi algoritmi che analizzeremo. Nello specifico, andremo a dividere il dataset in modo 'stratificato', al fine di avere la stessa distribuzione della variabile ' Retain' nei due dataset che otteniamo.

In [67]:
import itertools
def metrics(true, preds):
    """
    Function to calculate evaluation metrics 
    parameters: true values, predictions
    prints accuracy, recall, precision and f1 scores
    """

    accuracy = accuracy_score(true, preds)
    recall = recall_score(true, preds)
    precision = precision_score(true, preds)
    f1score = f1_score(true, preds)
    f2score= fbeta_score(true, preds,  beta=2, pos_label=1)
    print ('accuracy: {}, recall: {}, precision: {}, f1-score: {}, f2-score: {}'.format(accuracy, recall, precision, f1score, f2score))
    
    
def plot_confusion_matrix(cm,
                          target_names,
                          title='Matrice di confusione',
                          cmap=None,
                          normalize=True):
   

    accuracy = np.trace(cm) / np.sum(cm).astype('float')
    misclass = 1 - accuracy

    if cmap is None:
        cmap = plt.get_cmap('Blues')

    plt.figure(figsize=(8, 6))
    plt.imshow(cm, interpolation='nearest', cmap=cmap)#, aspect='auto')
    plt.title(title, fontsize=20)
    plt.colorbar()

    if target_names is not None:
        tick_marks = np.arange(len(target_names))
        plt.xticks(tick_marks, target_names, rotation=45, fontsize=15)
        plt.yticks(tick_marks, target_names, fontsize=15)

    if normalize:
        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]


    thresh = cm.max() / 1.5 if normalize else cm.max() / 2
    for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
        if normalize:
            plt.text(j, i, "{:0.4f}".format(cm[i, j]),
                     horizontalalignment="center",
                     color="white" if cm[i, j] > thresh else "black")
        else:
            plt.text(j, i, "{:,}".format(cm[i, j]),
                     horizontalalignment="center",
                     color="white" if cm[i, j] > thresh else "black")


    plt.tight_layout()
    plt.ylabel('True label', fontsize=15)
    plt.xlabel('Predicted label\n misclass={:0.4f}'.format(misclass), fontsize=15)
    plt.show()
  
In [69]:
X= data_clas.drop(['Retain','Profit0', 'District_first','AgeZero_0fascia','IncomeZero_0fascia'], axis=1)
y= data_clas['Retain']
x_train, x_test, y_train, y_test = train_test_split(X, y,
                                                    stratify=y, 
                                                    test_size=0.3, random_state=2)

X_TRAIN= x_train
Y_TRAIN=y_train

Regressione Logistica

Per la regressione logistica, sono stati seguiti due diversi approcci:

Primo approccio:

  1. Si ricercano i pesi della regressione logistica tramite cross validation ( il valore ottimale dei pesi è stato scelto tramite f1-score, che bilancia precisione and recall).
  2. Si performa la regressione logistica su tutto il dataset train con i pesi ottenuti.
  3. Si calcola la matrice di confusione sul dataset test.

Secondo approccio

  1. Si performa una regressione logistica pesata, in cui i pesi sono fissati tramite il comando class_weight='balanced'
  2. Si sceglie il valore di soglia della probabilità ( oltre la quale assegnare 1 ) in due possibili modi, cioè viene scelto :
    • il valore di probabilità che massimizza il G-mean ( una metrica che bilancia recall and specificità) nel seguente modo: $$G-mean= \sqrt{recall \cdot\ specificità}.$$ Se i positivi sono classicati erroneamente (cioè sono predetti come negativi), il valore G-mean è basso anche se i negativi sono classificati correttamente.
    • il valore di probabilità che massimizza F2-Measure ( una metrica che pone maggiore attenzione a minimizzare i falsi negativi ) , perchè è definita come $$ F2-Measure = \frac{5 \cdot Precision \cdot Recall}{4 \cdot Precision + Recall}$$.
  3. si calcola la matrice di confusione sul dataset test.

Regressione Logistica- I approccio

In [70]:
w = [ {0:1.0,1:1.0}, 
     {0:1.0,1:0.1}, {0:10,1:0.1}, {0:100,1:0.1},{0: 1.0, 1: 5} , {0: 1.0, 1: 9},
     {0:10,1:0.01}, {0:1.0,1:0.01}, {0:1.0,1:0.001}, {0:1.0,1:0.005}, 
     {0:1.0,1:10},{0:1.0,1:50}, {0:1.0,1:99}, {0:1.0,1:100},
     {0:1.0,1:200}, {0:1.0,1:1000}, {0:10,1:1000},{0:100,1:1000} ]
hyperparam_grid = {"class_weight": w}
lg3 = LogisticRegression(random_state=13)
grid = GridSearchCV(lg3,hyperparam_grid,scoring="f1", cv=100, n_jobs=-1, refit=True)
grid.fit(x_train,y_train)
print(f'Miglior valore di f1-score: {grid.best_score_} con parametri: {grid.best_params_}')
Miglior valore di f1-score: 0.4976939081395928 con parametri: {'class_weight': {0: 1.0, 1: 5}}

Performiamo la regressione logistica, con i pesi trovati dalla cross-validation e calcoliamo la matrice di confusione e alcune metriche sul dataset test.

In [71]:
lg3 = LogisticRegression(random_state=13, class_weight={0: 1.0, 1: 5})

#fit it
lg3.fit(x_train,y_train)
y_pred = lg3.predict(x_test)
plt.rcParams["axes.grid"] = False
lr_probs = lg3.predict_proba(x_test)
lr_probs = lr_probs[:, 1]

cm=confusion_matrix(y_test, y_pred)
target_names=['Clienti fedeli','Clienti persi']

plot_confusion_matrix(cm,
                          target_names,
                          title='Confusion matrix',
                          cmap=None,
                          normalize=False)
a=metrics(y_test, y_pred)
a

scores = pd.DataFrame(columns=['A'], index=['Regressione Logistica I','Regressione Logistica II','Regressione Logistica III','SVM','KNN'])

lr_probs = lg3.predict_proba(x_test)
lr_probs = lr_probs[:, 1]

forest_auc =  roc_auc_score(y_test, lr_probs)
forest_fpr, forest_tpr, _ =  roc_curve(y_test, lr_probs)

scores['A']['Regressione Logistica I'] =  forest_fpr, forest_tpr, forest_auc
accuracy: 0.7788431145295543, recall: 0.6838422391857506, precision: 0.4015689204333209, f1-score: 0.506001412096964, f2-score: 0.5995538204127162
In [73]:
metriche = pd.DataFrame(columns=['Recall', 'Misclass'], index=['Regressione Logistica I','Regressione Logistica II','Regressione Logistica III','SVM','KNN'])
misclass=1- accuracy_score(y_test, y_pred)
recall = recall_score(y_test, y_pred)
metriche['Recall']['Regressione Logistica I']= recall
metriche['Misclass']['Regressione Logistica I']= misclass

Regressione Logistica- II approccio

(G- mean)

In [74]:
from matplotlib import pyplot
from numpy import sqrt, argmax
model = LogisticRegression(solver='lbfgs', random_state=1, class_weight='balanced')
model.fit(x_train,y_train)
pyplot.rcParams["axes.grid"] = True

# predict probabilities
yhat = model.predict_proba(x_train)
# keep probabilities for the positive outcome only
yhat = yhat[:, 1]
# calculate roc curves
fpr, tpr, thresholds = roc_curve(y_train, yhat)
# calculate the g-mean for each threshold
gmeans = sqrt(tpr * (1-fpr))
# locate the index of the largest g-mean
ix = argmax(gmeans)
# plot the roc curve for the model
pyplot.plot([0,1], [0,1], linestyle='--', label='No-skill')
pyplot.plot(fpr, tpr,  label='Logistic')
pyplot.scatter(fpr[ix], tpr[ix], marker='o', color='black', label='Best')
# axis labels

pyplot.title('ROC Curve', fontsize=20)
pyplot.xlabel('False Positive Rate', fontsize=20)
pyplot.ylabel('True Positive Rate',fontsize=20)
pyplot.legend(prop={'size': 20})
# show the plot
pyplot.show()
print('Il valore di probabilità che massimizza G-Mean  è %f, con G-Mean=%.3f' % (thresholds[ix], gmeans[ix]))

model.fit(x_train,y_train)


y_pred = model.predict_proba(x_test)[:,1]
threshold= thresholds[ix]
y_pred = np.where(y_pred > threshold, 1, 0)
# Find prediction to the dataframe applying threshold


cm=confusion_matrix(y_test, y_pred)
target_names=['Clienti fedeli','Clienti persi']
plt.rcParams["axes.grid"] = False

plot_confusion_matrix(cm,
                          target_names,
                          title='Confusion matrix',
                          cmap=None,
                          normalize=False)
b=metrics(y_test, y_pred)
b

lr_probs = model.predict_proba(x_test)
lr_probs = lr_probs[:, 1]

forest_auc =  roc_auc_score(y_test, lr_probs)
forest_fpr, forest_tpr, _ =  roc_curve(y_test, lr_probs)

scores['A']['Regressione Logistica II'] =  forest_fpr, forest_tpr, forest_auc
misclass=1- accuracy_score(y_test, y_pred)
recall = recall_score(y_test, y_pred)
metriche['Recall']['Regressione Logistica II']= recall
metriche['Misclass']['Regressione Logistica II']= misclass
Il valore di probabilità che massimizza G-Mean  è 0.595556, con G-Mean=0.732
accuracy: 0.7947529238225688, recall: 0.6539440203562341, precision: 0.42269736842105265, f1-score: 0.5134865134865135, f2-score: 0.5894495412844037

(F2-Score)

Proviamo a cambiare soglia di probabilità.

In [75]:
model = LogisticRegression(solver='lbfgs', random_state=1, class_weight='balanced')
model.fit(x_train,y_train)
pyplot.rcParams["axes.grid"] = True

# predict probabilities
yhat = model.predict_proba(x_train)
# keep probabilities for the positive outcome only
yhat = yhat[:, 1]
# calculate pr-curve
precision, recall, thresholds = precision_recall_curve(y_train, yhat)
# convert to f score
fscore = (5 * precision * recall) / (4*precision + recall+ 0.000000000000000000001) #per non avere zero al denominatore
# locate the index of the largest f score
ix = argmax(fscore)
# plot the roc curve for the model
no_skill = len(y_train[y_train==1]) / len(y_train)

pyplot.plot([0,1], [no_skill,no_skill], linestyle='--', label='No-skill')
pyplot.plot(recall, precision, label='Regressione Logistica')
pyplot.scatter(recall[ix], precision[ix], marker='o', color='black', label='Best')
# axis labels
pyplot.title('Curva Precisione vs Recall', fontsize=20)
pyplot.xlabel('Recall', fontsize=20)
pyplot.ylabel('Precision', fontsize=20)
pyplot.legend(prop={'size': 20})
# show the plot
pyplot.show()
print('Il valore  di probabilità che massimizza F2-score è %f, con F2-Score=%.3f' % (thresholds[ix], fscore[ix]))


model.fit(x_train,y_train)

y_pred = model.predict_proba(x_test)[:,1]
threshold= thresholds[ix]
y_pred = np.where(y_pred > threshold, 1, 0)
# Find prediction to the dataframe applying threshold


cm=confusion_matrix(y_test, y_pred)
target_names=['Clienti fedeli','Clienti persi']
plt.rcParams["axes.grid"] = False

plot_confusion_matrix(cm,
                          target_names,
                          title='Confusion matrix',
                          cmap=None,
                          normalize=False)
c=metrics(y_test, y_pred)
c
lr_probs = model.predict_proba(x_test)
lr_probs = lr_probs[:, 1]

forest_auc =  roc_auc_score(y_test, lr_probs)
forest_fpr, forest_tpr, _ =  roc_curve(y_test, lr_probs)

scores['A']['Regressione Logistica III'] =  forest_fpr, forest_tpr, forest_auc
misclass=1- accuracy_score(y_test, y_pred)
recall = recall_score(y_test, y_pred)
metriche['Recall']['Regressione Logistica III']= recall
metriche['Misclass']['Regressione Logistica III']= misclass
Il valore  di probabilità che massimizza F2-score è 0.406392, con F2-Score=0.597
accuracy: 0.7428089769255084, recall: 0.7245547073791349, precision: 0.3619319987289482, f1-score: 0.48272939182030095, f2-score: 0.6036036036036037

SVM

Per i successivi algoritmi, normalizziamo il dataset, in quanto sono algoritmi che si basano su concetti di distanza.

In [76]:
scaler = StandardScaler()  
a=  np.asarray( x_train.columns.tolist())
b=  np.asarray( x_test.columns.tolist())
X_TRAIN = scaler.fit_transform(X_TRAIN)
X_TRAIN=pd.DataFrame(X_TRAIN , columns=a )
x_test = scaler.fit_transform(x_test)
x_test=pd.DataFrame(x_test , columns=b)

Performiamo la cross-validation (K=5) per cercare i tuning parameters degli algoritmi.

In particolare, applichiamo una versione 'modificata' della cross-validation:

  • in ogni iterazione, è esclusa una parte del dataset per la validazione
  • le k-1 folds del training vengono bilanciate tramite undersampling ( RandomUnderSampler), mentre la fold Validation non viene modificata.
  • Ripetiamo K volte

In questo modo, i dati in fase di validazione sono rappresentativi del dataset reale.

In sintesi, il modello è 'allenato' su un dataset bilanciato, mentre la validazione avviene su un dataset non bilanciato.

Dopo aver scelto il valore ottimale del parametro, si fitta il modello su un dataset bilanciato. Quindi si calcola la matrice di confusione sul dataset test non bilanciato.

Il valore ottimo del tuning parametro è scelto sulla base del valore di F2-score. In particolare:

  • per SVM, cerchiamo il valore ottimo di C, che è un parametro di regolarizzazione. Come possibili valori di C, proviamo C=[0.1, 1, 10 , 100, 1000];
  • per KNN, cerchiamo il valore ottimo di K, che è il numero di vicini che vengono considerati dall'algoritmo. Come possibili valori di K, proviamo K da 1 a 20.
In [77]:
from imblearn.under_sampling import RandomUnderSampler
spl=5
random_seed=1
kf = StratifiedKFold(n_splits=spl, random_state=random_seed)
C_range = [0.1, 1, 10, 100, 1000]
griglia= 5
cross_val_fbeta_score_lst= np.zeros(shape=(len(C_range),spl))
rus = RandomUnderSampler(random_state=0, replacement=True)

I boxplot che seguono mostrano le statistiche su F2-score ottenute durante la cross-validation per ogni valore di C. Come vediamo, il migliore valore di C è 10.

In [78]:
#svm
i=-1
for train_index_ls, validation_index_ls in kf.split(X_TRAIN, Y_TRAIN):
    train, validation =X_TRAIN.iloc[train_index_ls], X_TRAIN.iloc[validation_index_ls]
    target_train, target_val =  Y_TRAIN.iloc[train_index_ls],  Y_TRAIN.iloc[validation_index_ls]
    rus = RandomUnderSampler(random_state=0, replacement=True)
    X_train_res, y_train_res = rus.fit_resample(train, target_train)
    i=i+1
   
    for C in range( len( C_range)) :
        model = SVC(C=C_range[C])
        model.fit(X_train_res, y_train_res)
        validation_preds= model.predict(validation)
        cross_val_fbeta_score_lst[C][i]=fbeta_score(target_val, validation_preds,  beta=2, pos_label=1)
In [79]:
df = pd.DataFrame(cross_val_fbeta_score_lst.transpose(),
                        columns= ['C=0.1','C=1', 'C=10', 'C=100', 'C=1000'])
    #boxplot = df.boxplot( rot=45, fontsize=15)
#df.plot.box(return_type='axes', figsize=(12,8))

boxplot = df.boxplot(column=['C=0.1','C=1', 'C=10', 'C=100', 'C=1000'], fontsize=15, rot=45, figsize=(12,8))

   # print ('Cross validated beta score: {}'.format(cross_val_fbeta_score_lst.mean(axis=1)))
vettore_medie=cross_val_fbeta_score_lst.mean(axis=1)
m=np.argmax(vettore_medie)
#cross_val_fbeta_score_lst

Quindi fittiamo il modello con C=1 su il dataset train, dopo averlo bilanciato.

In [80]:
x_train_rus, y_train_rus = rus.fit_resample(X_TRAIN, Y_TRAIN)
model = SVC(C=1, probability=True)
model.fit(x_train_rus, y_train_rus)
Out[80]:
SVC(C=1, cache_size=200, class_weight=None, coef0=0.0,
    decision_function_shape='ovr', degree=3, gamma='auto_deprecated',
    kernel='rbf', max_iter=-1, probability=True, random_state=None,
    shrinking=True, tol=0.001, verbose=False)

Riportiamo quindi le metriche e la matrice di confusione sul dataset test.

In [81]:
validation_preds= model.predict(x_test)


cm=confusion_matrix(y_test, validation_preds)
plt.rcParams["axes.grid"] = False

plot_confusion_matrix(cm,
                          target_names=['Clienti fedeli','Clienti persi'],
                          title='Confusion matrix',
                          cmap=None,
                          normalize=False)

d=metrics(y_test, validation_preds)
d


lr_probs = model.predict_proba(x_test)
lr_probs = lr_probs[:, 1]
forest_auc =  roc_auc_score(y_test, lr_probs)
forest_fpr, forest_tpr, _ =  roc_curve(y_test, lr_probs)

scores['A']['SVM'] =  forest_fpr, forest_tpr, forest_auc
misclass=1- accuracy_score(y_test, validation_preds)
recall = recall_score(y_test, validation_preds)
metriche['Recall']['SVM']= recall
metriche['Misclass']['SVM']= misclass
accuracy: 0.7866399747128859, recall: 0.6692111959287532, precision: 0.4114196323816973, f1-score: 0.5095664809881326, f2-score: 0.5946862634256641

KNN

Ricerchiamo il valore ottimo di K mediante cross-validation.

In [82]:
random_seed=1
spl=5
griglia=20
kf = StratifiedKFold(n_splits=spl, random_state=random_seed)
cross_val_fbeta_score_lst= np.zeros(shape=(griglia,spl))


#KNN
    # Train the model using the training sets
i=-1
for train_index_ls, validation_index_ls in kf.split(X_TRAIN, Y_TRAIN):
        # keeping validation set apart and oversampling in each iteration using smote 
    train, validation =X_TRAIN.iloc[train_index_ls], X_TRAIN.iloc[validation_index_ls]
    target_train, target_val =  Y_TRAIN.iloc[train_index_ls],  Y_TRAIN.iloc[validation_index_ls]
    rus = RandomUnderSampler(random_state=0, replacement=True)
    X_train_res, y_train_res = rus.fit_resample(X_TRAIN, Y_TRAIN)
    i=i+1
    for k in range(1,griglia+1): #k=1 a k=griglia
        model = KNeighborsClassifier(n_neighbors=k)
        model.fit(X_train_res, y_train_res)
        validation_preds= model.predict(validation)
        cross_val_fbeta_score_lst[k-1][i]=fbeta_score(target_val, validation_preds,  beta=2, pos_label=1)
        #cross_val_fbeta_score_lst[k-1][i]=recall_score(target_val, validation_preds, pos_label=1)


#print ('Cross validated beta score: {}'.format(cross_val_fbeta_score_lst.mean(axis=1)))
vettore_medie=cross_val_fbeta_score_lst.mean(axis=1)
train_mean = np.mean(cross_val_fbeta_score_lst, axis=1)
train_std = np.std(cross_val_fbeta_score_lst, axis=1)
# Draw lines
plt.plot(range(1,griglia+1), train_mean, '--', color="#111111",  label="Training score")

# Draw bands
plt.fill_between(range(1,griglia+1), train_mean - train_std, train_mean + train_std, color="#DDDDDD")
# Create plot
plt.title("cross validation Curve", fontsize=20)
plt.xlabel("K", fontsize=20), plt.ylabel("2-SCORE",  fontsize=20), plt.legend(loc="best", prop={'size': 20})
plt.tight_layout()

m=np.argmax(vettore_medie)
k_ottimo=m+1 #nella posizione zero del vettore c'è k=1
plt.plot(k_ottimo,max(vettore_medie),'ro') 
plt.xticks( range(1,griglia+1)  ,range(1,griglia+1), fontsize=20)
plt.yticks(fontsize=20)
plt.show()
print('Il valore ottimale di K è 1.')
Il valore ottimale di K è 1.

Alleniamo il modello sul dataset train bilanciato e riportiamo le metriche e la matrice di confusione sul dataset test.

In [83]:
x_train_rus, y_train_rus = rus.fit_resample(X_TRAIN, Y_TRAIN)

model = KNeighborsClassifier(n_neighbors=1)
model.fit(x_train_rus, y_train_rus)

# F0.5-measure puts more attention on minimizing false positives t
#positivo= 1 
#voglio minimizzare chi è zero ma dice che è 1 , cioè falsi negativi
validation_preds= model.predict(x_test)




cm=confusion_matrix(y_test, validation_preds)
plt.rcParams["axes.grid"] = False

plot_confusion_matrix(cm,
                          target_names=['Clienti fedeli','Clienti persi'],
                          title='Confusion matrix',
                          cmap=None,
                          normalize=False)

e= metrics(y_test, validation_preds)
e
lr_probs = model.predict_proba(x_test)
lr_probs = lr_probs[:, 1]
forest_auc =  roc_auc_score(y_test, lr_probs)
forest_fpr, forest_tpr, _ =  roc_curve(y_test, lr_probs)

scores['A']['KNN'] = forest_fpr, forest_tpr, forest_auc  
misclass=1- accuracy_score(y_test, validation_preds)
recall = recall_score(y_test, validation_preds)
metriche['Recall']['KNN']= recall
metriche['Misclass']['KNN']= misclass
accuracy: 0.6152144136550416, recall: 0.6552162849872774, precision: 0.24879227053140096, f1-score: 0.36064425770308123, f2-score: 0.4938626774069812

Conclusioni

Abbiamo visto 5 diversi modelli di classificazione.

Il peggior modello è il KNN. I modelli di regressione logistica risultano essere, in generale, i migliori.

Per riassumere il lavoro di classificazione mostriamo in seguito il valore di misclass, ( 1- accuracy, cioè la percentuale di osservazioni mal classificate) , il valore di Recall ( che vale 1 se non ci sono falsi negativi) e la ROC Curve.

In [84]:
fig = go.Figure()

fig = fig.add_trace(go.Scatter(x=scores['A']['Regressione Logistica I'][0], y=scores['A']['Regressione Logistica I'][1], mode='lines',
                                       line_shape='linear', name='Regressione Logistica I'))    
fig = fig.add_trace(go.Scatter(x=scores['A']['Regressione Logistica II'][0], y=scores['A']['Regressione Logistica II'][1], mode='lines',
                                       line_shape='linear', name='Regressione Logistica II'))
fig = fig.add_trace(go.Scatter(x=scores['A']['Regressione Logistica III'][0], y=scores['A']['Regressione Logistica III'][1], mode='lines',
                                       line_shape='linear', name='Regressione Logistica III'))
fig = fig.add_trace(go.Scatter(x=scores['A']['SVM'][0], y=scores['A']['SVM'][1], mode='lines',
                                       line_shape='linear', name='SVM'))
fig = fig.add_trace(go.Scatter(x=scores['A']['KNN'][0], y=scores['A']['KNN'][1], mode='lines',
                                       line_shape='linear', name='KNN'))

fig.update_xaxes(range=[-0.05, 1.05], title_text='False Positive Rate')
fig.update_yaxes(title_text='True Positive Rate')
fig.update_layout(title='ROC Curve')
fig.show()
In [85]:
metriche
Out[85]:
Recall Misclass
Regressione Logistica I 0.683842 0.221157
Regressione Logistica II 0.653944 0.205247
Regressione Logistica III 0.724555 0.257191
SVM 0.669211 0.21336
KNN 0.655216 0.384786